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ABSTRACT 

We calculate the probability that a Milky-Way-like halo in the standard cosmological model has the observed 
number of Magellanic Clouds (MCs). The statistics of the number of MCs in the ACDM model are in good 
agreement with observations of a large sample of SDSS galaxies. Under the sub-halo abundance matching 
assumption of a relationship with small scatter between galaxy r-band luminosities an d halo internal ve locities 
Vmax, we make detailed comparisons to similar measurements using SDSS DR7 data bv lLiu et "ail ( 120101) . Mod- 
els and observational data give very similar probabilities for having zero, one, and two MC-like satellites. In 
both cases, Milky Way-luminosity hosts have just a ~ 10% chance of hosting two satellites similar to the Mag- 
ellanic Clouds. In addition, we present a prediction for the probability for a host galaxy to have N sats satellite 
galaxies as a function of the magnitudes of both the host and satellite. This probability and its scaling with host 
properties is significantly different from that of mass-selected objects because of scatter in the mass-luminosity 
relation and because of variations in the star formation efficiency with halo mass. 

Subject headings: galaxies: dwarf, Magellanic Clouds, evolution — dark matter 



1. INTRODUCTION 

Understanding the satellite population of our own Milky 
Way galaxy is one of the outstanding problems in galaxy evo- 
lution studies today. Because these objects have weak grav- 
itational potential wells, it is quite likely that different phys- 
ical processes played a driving role in determining how such 
galaxies formed and evolved than for brighter, more massive 
objects. One manifestation of this is the well-known "miss- 
ing satellites problem," which refers to the fact that the num- 
ber of satellite galaxies predicted from simulations appears 
to be much high er than the number actually observed a round 
the Milky Way dKlypin et al.lll999b iMoore et al.lll999l) . One 
common way to address this problem has been to run simu- 
lations of individual objects using extremely high mass reso- 
lution. To date, roughly 10 of these ultra-high resolution N- 
body simulations have been run, which typically resolve indi- 
vidual Milky Way mass halos wit h upwards of 10 9 particl es 
dMadau et al.H2008t LSpringel et al.1l2008t LStadel et al.ll2009h . 

The Magellanic Clouds (MCs) are a pair of well studied 
satellites of the Milky Way (MW) in the Southern Hemi- 
sphere with magnitudes My = -18.5 and -17.1 for the Large 
and Small Magellanic Clouds (LMC and SMC). While their 
large angular size makes detailed mass measurements some- 
what difficult, observations constrain th eir maximum circular 
velocities to be arou nd v mj i , ~ 60 km s~ 1 dvan den Berghl2000t 
Ivan der Marel et al.ll2002f1stanimirovic et al.l 120041) . Unlike 
the missing satellite problem, one common feature of these 
ultra-high resolution simulations is that they typically under- 
predict the number of massive sa tellites compared to the MCs 
in the MW (see, e.g., Figure 5 in lMadau et al.ll2008l) . Indeed, 
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most of these high resolution halos contain no objects of sim- 
ilar mass to either of the two MCs. While this should hardly 
be taken as a serious problem since the handful of high resolu- 
tion simulated halos does not constitute a statistical sample, it 
is important to understand how typical or atypical the Milky 
Way is, and whether or not not there is an "extra satellites" 
problem at high satellite masses within the modern Lambda 
Cold Dark Matter (ACDM) paradigm. 

Because the MCs are relatively nearby and have been stud- 
ied with high-precision instruments such as HST for some 
time, detailed 3-dime nsional measurements of their veloci- 
ties have been made dKallivayalil et al.l 2006b a). Some of 
these studies have indicated that the Magel lanic Clouds may 
be on their first or bit around the Milky Way (iBesla et al.l2007t 
iBusha et al.l2010t) : this may make such an extra satellite prob- 
lem less worrisome if the presence of the Magellanic Clouds 
is a transient event. Regardless of the dynamical state of 
the Clouds, it is an important test of galaxy formation to un- 
derstand the likelihood for MW-like systems to host massive 
satellite galaxies. 

Recent developments have afforded us the possibility to 
address this problem in more detail with both observations 
and theoretical models. From the observational side, wide 
area surveys such as the Sloan Digit al Sky Survey (SDSS 
I York et al.l20"00l:lAbazajian et al.l2009l) have given us the abil- 
ity to probe galaxy c ontent for a large number of Milky Way- 
magnitude obj ects (IChen et al.l l2006t iJames & Ivorvl 120101 : 
iLiu et aUHoTol) . 

The main theoretical effort has been in populating dark mat- 
ter halos with galaxies u sing semi-analy t ic me thods and hy- 
drodynamic simulations. iKoposov et al.1 (120091) used a num- 
ber of toy models to add galaxy properties to dark matter 
halos generated using a combination of Press-Schechter the- 
ory with semi-analy tic models for tracking subhalo orbits 
dZentner et al.ll2005l) . They found that it was very difficult to 
model objects as bright as Magellanic clouds without allow- 
ing for an extremely high star formation efficiency. This find- 
ing extended the underprediction of high-mass subhalos to an 
underprediction of luminous satellites, implying that the MW- 
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MCs system is unusual. Similarly, Okamoto et al.l (1201 Oh ex- 
plored a range of feedback models to add galaxies to som e of 
the high resolution Aquarius halos dSpringel et al. 2008|). It 
was again difficult to readily reproduce halos with luminosi- 
ties as bright as the MCs. Having looked in detail at a handful 
of simulated objects, however, this work indicates that there 
may be significant halo-to-halo scatter in the number of such 
massive objects. The idea of intrinsic scatter in the subhalo 
population was expanded on in Ishi vama et al.ld2009al) . While 
this work did not concentrate specifically on MC-like subha- 
los, they did consider the range in number of subhalos with 
Vmax,sub/vmaxjwsr > 0.1. This work showed an extremely large 
variation (20-60) in the number of such massive subhalos a 
galaxy-sized halos would host. 

In contrast to the se mi-analytic modeling just discussed, 
iLibeskind et al.l (120071) used hydrodynamic simulations to 
model the luminosity functions for the satellites around MW- 
like host halos. Their simulation identified 9 MW-like central 
galaxies and found that they, on average, have 1.6 satellites 
brighter than My = -16 and a third of them had satellites with 
luminosities comparable to the LMC. 

The recent Millennium - II and Bol s hoi simulations 
dBoylan-Kolchin et al.l 120101: iKlypin et al.l 120101) have, for 
the first time, allowed us to probe cosmological volumes to 
understand the ACDM predictions for the satellite popula- 
tions of MW-like halos. The p roperties of these simulations 
are summarized in Table 1. Bovlan-Kolc hin et~ail ( 120101) 
(hereafter BK10) quantified the likelihood for 10 12 M Q 
halos to host massive satellite galaxies in the Millennium-II 
simulation, finding that subhalos similar to the MCs are 
quite rare. In this paper, we expand on this work by making 
similar measurements for the Bolshoi simulation, which used 
WMAP7 cosmological parameters, and using an abundance 
matching technique to make detailed comparisons between 
the Bo lshoi predictions and the measurements from lLiu et al.l 
( 120101) . Our goal is to understand just how well ACDM 
reproduces the statistical properties of bright satellites. 

Note that this is the revers e question from the one that was 
asked in a companion paper. iBusha et al] ( 120101) . That work 
assumed a satellite population and asked what the implica- 
tions were for the properties of the host halo, including its 
mass. Here, we assume a host halo mass and ask about the 
implications for the subhalo population. There is no reason 
for bo th questions to give the same answer: while lBusha et alj 
(120101) showed that a halo which hosts two MC-like satellites 
most likely has a mass near 1.2 x 10 12 M Q , there is no reason 
to assume that a typical 1.2 x 10 12 M© halo will have the MCs 
as satellites. 

We begin by giving an overview of the Bolshoi simulation 
in §|2] and then investigate the properties of massive dark mat- 
ter satellites around dark matter host halos in §[3] focusing 
on the mass ranges for hosts and satellites that are most rel- 
evant to the MW system. The analysis here is similar to that 
of BK10. In §H1 we assign galaxy luminosities to our suite of 
dark matter halos and extend the results for a sample with sim- 
ilar selection cuts as for observations. In this way, we are able 
to make detaile d comparisons to the observational work of 
dLiu et al.ll2010l hereafter LI 0) concerning the s atell ite popu- 
lation around MW-magnitude galaxies. Section l4~4l gives the 
results of this analysis — see especially Figure [8] Finally, 
in §|5] we expand this study to include the satellite popula- 
tion of a more general distribution of hosts, and §|6] summa- 
rizes our conclusions. Throughout this paper, we adopt the 



convention h = 0.7 (the value that was used in the Bolshoi 
simulation) when reporting values from either simulations or 
observations. 

2. SIMULATIONS 

We u se the dark matter halos identified in the Bolsho i sim- 
ulation dKlvpin et al.ll2010tlTruiillo-Gomez et alj|2010l) . This 
simulation modeled a 250 /r'Mpc comoving box using cos- 
mological parameters similar to those derived by WMAP7 
dKomatsu et al.l |2010|) : Q m = 0.27, ft A = 0.73, er 8 = 0.82, 
n = 0.95, and h = 0.7. The simulation volume contains 2048 3 
particles, each with a mass of 1.35 x 10 8 /r 'M,^ and was run 
using the ART code dKravtsov et al.l 1 19971) . 180 snapshots 
from the simulation were saved and analyzed. One of the 
unique aspects of this simulation is the high level of spatial 
resolution employed, allowing objects to be resolved down to 
a physical scale of 1 /i _1 kpc. This gives us excellent ability to 
track halos as they merge with and are disrupted by larger ob- 
jects, allowing us to track them even as they pass near the core 
of the host halo. A summary of these simulation parameters 
is presented in Table Q] Because we discuss our work in the 
context of the BK10 results, we also pr esent the same param- 
eters for the Millennium II simulation dBoylan-Kolchin et al.l 
120091) . on which the BK10 results were based. 

Hal os and subhalos were ident ified using the BDM algo- 
rithm dKlvpin & Holtzmanlll997l) . The algorithm identifies 
maxima in the density field and examines the neighboring re- 
gion to identify bound particles. In this way, it treats both 
halos and subhalos identically. Subhalos are just identified 
as objects living within the virial radius of a larger objects. 
Because of the high level of mass and spatial resolution, 
BDM results in a halo catalog that is complete down to a 
maximum circular velocity v max = 50 km s -1 , where v max = 

max fiy/ GM(< r)/rj . This corresponds to a virial mass of 

roughly 10 I0 /r I M Q . When BDM halos are identified, the ids 
of their 50 most bound particles are also stored to assist in 
producing merger trees. 

As discussed in Section l4~2l in order to assign galaxy lumi- 
nosities to dark matter halos, we need to track the histories of 
dark matter substructures. This is done using merger trees cre- 
ated from the 180 snapshots of the Bolshoi simulation. The 
de tailed algorithm for cre ating the merger trees is described 
in iBehroozi et al.l d2010al) . Briefly, the algorithm works by 
first linking halos across time steps by tracking the 50 most 
bound particles of each halo. Some halos will not have any 
of their 50 most bound particles identified at a later timestep 
(e.g., very massive objects in which the 50 most bound parti- 
cles change rapidly through stochastic processes), while some 
will have their particles distributed to multiple halos. The al- 
gorithm corrects for this by running a simple N-body calcula- 
tion on the locations and masses of all halos in the simulation 
to predict where each halo should wind up at the next time 
step. Using this information, it is possible to link halos across 
multiple time steps with very high accuracy. 

3. SATELLITE STATISTICS FOR DARK MATTER SUBHALOS 

We begin by investigating the properties of dark matter 
satellites around dark matter hosts in the Bolshoi simulation, 
focusing on the mass range for host s and satellites that is rel- 
evant to the MW syste m. In 33.11 we consider trends with 
host halo mass, and in 33.21 we investigate trends with halo 
environm ent. W e then consider the distribution of the satellite 
number ( 33.3b . Finally, we extend this analysis to include a 
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TABLE 1 

Comparison of the simulation parameters from Bolshoi (used 
for this work) with those of the millennium ii simulation 

(USEDINBK10) 



Bolshoi 



MSII 



Box Size [Mpc] 
Np 
M P [M @ ] 
h 

s hi 



357 
2048 3 
1.9 x 10 8 
0.7 
0.27 
0.82 



143 
2160 3 
9.8 x 10 6 
0.73 
0.25 
0.9 



Force Resolution [kpc] 1 .4 (proper) 1 .4 (Plummer-equivalent softening) 

more observationally relevant selection based on galaxy lumi- 
nosities in $U 

Because the mass resolution of the Bolshoi simula- 
tion creates a halo catalog complete down to 50 km s" 1 , 
roughly equivalent to the lower bound of the mass of 
the MCs (Ivan der Marel et al.ll2002t IStanimirovic et alJl2004t 
lHarris & Zaritsky||2006l) . we begin by measuring the proba- 
bility distribution for halos hosting N su bs subhalos with v max 
larger than a given value. This measurement is made by 
identifying the virial radius of all distinct (non-subhalo) MW 
mass halos and counting the number of objects internal to 
their virial radii. We take the virial mass of the Milky Way 
to be log(M vir /M ) = 12.08 ±0.1 2 (Mw = 1.2 x 10 12 M Q ) 



the mass measured in iBusha et"afl (12010 ), which is consis- 
tent with most results in the literature (iBattaglia et all 12005b 
iDehnen et al1l2006t ISmith et alj|2007t IXue et al.ll2008l This 
results in 36,000 MW analogs found in the Bolshoi volume. 
We present the resulting probability distribution of satellite 
counts as the colored points in the left panel of FigureQ] Here, 
the different colors represent different thresholds for satellite 
Vmax: red, green, and blue represent all satellites with v max > 
50, 60, and 70 km s" 1 . Error bars were calculated using the 
bootstrap method and represent 95% confidence intervals. 

From Figure Q] we can immediately see that the likelihood 
of hosting multiple satellites is somewhat low. For satel- 
lites with v max greater than 50 km s" 1 , roughly 30% of sim- 
ulated hosts contain one or more subhalos and just 8% have 
two or more. These numbers are in excellent agreement with 
the work of BK10, who performed a similar analysis on the 
Millennium II simulation, which was run using a different N- 
body code in a WMAP1 cosmology and with a very different 
subhalo identification algorithm. The probability of hosting 
multiple satellites drops precipitously with increasing satel- 
lite mass; e.g., the probability of a halo hosting two subhalos 
larger than v max = 70km s" 1 is just 1 %. 

For comparis on, the black poin ts in FigureQ]represent mea- 
surements from iLiu et al.l (120101) . who measured the proba- 
bility distribution for finding MC-luminosity satellites within 
250 kpc around isolated MW-luminosity galaxies (the mean 
virial radius of our sample). We only plot these measure- 
ments out to N su bs = 3 because, for a 250 kpc aperture, uncer- 
tainties from their background subtraction become extremely 
large for higher (lower probability) values of N s „bs- Here, the 
agreement for satellites larger than 50 km s _1 is striking, in- 
dicating that CDM can reproduce the statistics of MCs. We 
must, however, be careful about over-interpreting this result. 
In particular, their selection criteria are very different from 
our mass selection. We will return to this in §4.4 where we 



make a comparison directly to similarly selected samples. 

By comparison with the Bolshoi simulation, the existence 
of two satellites with v max larger than 50km s -1 makes the 
MW almost a 2a outlier. Understanding whether anything 
other than random chance is responsible for putting the MW 
in this slightly rare regime motivates further investigation into 
which other properties may correlate with satellite abundance. 
In particular, we investigate the correlation of the number of 
subhalos with host mass and environment. 



3.1. Mass dependence 

The correlation between sub halo abundance an d 
mass has been we l l studied (IZentner & Bullock! 
Diemand et all 120041: iGao et al.l 120041: IZentner et al.l 



host 



2003 



2005 



Klypin et al These studies have all found that, when 



scaled in units of fi = M su b/Mi wst , larger halos contain more 
subhalos above a given p than do smaller ones. This is 
the expected result from the hierarchical merging picture of 
CDM, in which larger objects grow throu gh the merging of 
smaller objects and there fore form later dBlumenthal et al.l 
Il984t lLacey & Colell 19931) . This later formation has a two- 
fold impact: it reduces the amount of time available to tidally 
strip a subhalo, and lowers the host concentration, further 
reducing the ability of a host to tidall y strip its subhalos 
(IWechsler et al.ll2002b IBusha et al.ll2007l) . Motivated by this, 
we present the abundance of LMC and SMC-mass subhalos 
for the 27,000 Bolshoi halos with M vir = 2.6 ±0.9 x 1O 12 M . 
The results can be seen in FigureQ] Here, the left panel shows 
the results for selecting hosts with M vir = 0.8-1.7 x 10 12 M Q , 
while the right panel shows a mass range M wu = 1.7-3.4 x 
10 12 M Q . 

For the more massive M V1I ~ 2.6 x 10 12 M Q halos, the cu- 
mulative probability of hosting two or more 50 km s" 1 halos 
increases by a factor of 3 compared to the 1.2 x 10 12 M Q ob- 
jects (from 8% to 26%), while jumping by a factor of 6 for 
70km s" 1 subhalos (from 1% to nearly 6%). W e emphasize 
that th is is not a contradiction with the results of IBusha et all 
(120101) . which found that a halo with exactly two MC-like 
subhalos was likely to have mass 1.2 x 1O 12 M . While the 
fraction of more massive halos hosting two MC-like subha- 
los is larger, the steep slope of the mass function means that 
there are many more lowe r-mass halos . Add itionally, aside 
from just considering v max , IBusha et"ai1 (1201 Oh selected halos 
with satellites like the MCs in terms of location and three- 
dimensional velocity, properties not considered in the present 
analysis. 

Again, these results are in excellent agreement with the 
work of BK10. While such an agreement is generally 
expected from pure collisionless N-body simulations, such 
agreement is not as trivial as it may appear. While the present 
work uses results from the Bolshoi simulation, run using 
the adaptive refinement ART code and the BDM halo finder, 
BK10 used the Millennium-II simulation run with the TreePM 
code Gadg et-2 with substructur es identified using the Subfind 
algorithm (ISpringel et al.ll200ll) . Additionally, in both cases, 
we are concerned with subhalos close to the mass-limit for 
the simulations. Thus, the close agreement between these two 
works suggests that numerical effects in the simulations are 
well understood at this level. 

3.2. Environmental dependence 

Beyond simple host mass dependence, the size of the cos- 
mological volume modeled by the Bolshoi simulation allows 
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FIG. 1 . — The probability distribution for the abundance of MC-like satellites around MW-mass hosts in simulations (colored diamonds). The left panel shows 
the results for hosts in a mass bin of M vu = 1.2 ±0.3 X 10 i2 Mq; the right panel shows results for hosts in a mass bin of M vir = 2.6 ± 0.4 X 1O 12 M . In both 
cases, the red, green, and blue distributions correspond to the number of satellites with v max > 50km s , Vmax > 60km s~', and v max > 70km s _I , respectively. 
Dashed lines show th e best-fit negative binomial distribution, and dotted lines show the best-fit Poisson distribution. In both panels, the black points show the 
observational result of Liu et al. (2010) for the number of satellites within 250 kpc around isolated MW-luminous galaxies. 



us to perform further detailed studies on properties impact- 
ing the satellite abundance distribution. Here, we quantify the 
environmental dependence of the subhalo population by split- 
ting the sample on local density, defined by 



Ph.r-Ph 
Ph 



(1) 



where /?/, ,. is the Eulerian density of dark matter contained in 
halos larger than M vil =2x10' 'M (roughly 1000 particles in 
Bolshoi) defined in a sphere of size r. In the present analysis, 
we take r = 1.0/i _1 Mpc; we have performed similar analysis 
on a range of scales and find this to be the range where the 
galaxy halo distribution is the most sensitive to the presence 
of massive satellites. 

We compare the distribution P(N su i, s ) for halos in the top 
and bottom quartile distributions to the mean relation in Fig- 
ure I2 There is a clear systematic trend: halos living in over- 
dense regions (red points) are more likely to host massive sub- 
halos than those in underdense regions (blue points). This is 
likely due to the fact that halos in denser environments live in 
more biased regions, where smaller density perturbations are 
amplified, making it easier for smaller halos to form and ac- 
crete onto intermediate mass hosts. This makes the deviations 
from the mean relation in Figure |2] a potentially observable 
manifestation of assembly bias, which posits that the prop- 
erties of a given halo may be determined by properties other 
than mass, such as the halo environment. 

Selecting for objects in a dense environment increases by 
almost 25% the probability of having N su b s = 2 or more, from 
9% to just over 11%; this effect becomes stronger for larger 
values of N su bs ■ The Milky Way is somewhat overdense on this 
scale, having a massive companion, M31, as well as a num- 
ber of surrounding dwarf galaxies. Recent estimates put the 
total m ass of the local group at roughly 5 x 10 12 M Q dLi et al.l 
l2009l) . solidly in a dense environment on the ~ 1 Mpc scale. 
Even if we ignore the contribution to the local density from all 
neighboring objects except M31, the MW lies close to the top 
quartile of our measured local density distribution in Bolshoi 
on a ~ 1 Mpc scale. We can begin to directly quantify the im- 
pact of such a massive neighbor on the subhalo population by 
identifying MW-mass halos which have M31-like neighbors 
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FIG. 2. — The probability that a 1.2 X 1O 12 M0 halo hosts N su b s subhalos 
with v max > 50km s , split by local host density. The black points represent 
all halos, while the red and blue points represent the most overdense and 
underdense quartiles, respectively. Green circles represent MW-mass halos 
that were selected to have M31-mass neighbors. 



with M vir = (1 -3) x 1O 12 M at a distance of 700-900 kpc. 
These systems are represented by the green circles in Figure 
|2j which have a N su b s = 2 probability enhanced by la relative 
to the full sample. Thus, the local environment around the 
MW, which is dominated by the presence of M31, boosts the 
chance of seeing the Magellanic Clouds in the Milky Way. 

3.3. Probability distributions 

As an extension of the trends shown in Figure Q] it is use- 
ful to generalize these results to arbitrary host and subhalo 
masses by fitting the probability distribution. A Poisson dis- 
tribution seems to be the most natural assu mption, and has 
been found previously in the literature (e.g.. iKravtsov et al] 
120041) . Deviations from a Poisson distribution can be quanti- 
fied by the calculation of the second moment, oii, where 



Oil = 



(N(N-l)) 
(AO 2 ' 



(2) 
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For a Poisson distribution, oli = 1. While studies have found 
some evidence for a deviation from Poisso n, with 02 > 1 
dKravtsov et al.ll2004t IWetzel & Whitell2010l) . the distribution 
has been more fully quantified by BK10. This work found 
that the distribution could be much better modeled using a 
negative binomial distribution, 



P(N\r,p): 



T(N+r) 

r<»r(jv+ 1) 



P r a- P f 



where 



P = 



1 



l+(a 2 2 -l)(N) 



1 



(3) 



(4) 



although they note that a Poisson distribution is a reasonable 
fit when (N) is low, such as the case of the MCs. We plot 
fits to both a Poisson and a negative binomial distribution to 
our measurements as the dotted and dashed lines in Figure 
Q] We again see that, because (N) < 1, the Poisson distri- 
bution provides a reasonable fit, while the negative binomial 
distribution fit is excellent (at the expense of adding an addi- 
tional parameter). 6 Our results confirm those of BK10, in that 
the distribution becomes significantly more non-Poisson for 
larger values of (N). 

As a final comparison of the Bolshoi halos to previous work 
in the literature, we plot in Figure [3] the trends of (N) (top 
panel, similar information to Figure 15 in lKlvpin et aUl20lol) 
and «2 (bottom panel) with mass ratio v acC:SU b /v max ^ost for a 
range of host masses. Here, v acc , s „b is the maximum circu- 
lar velocity, v max , that the subhalo had at the epoch of ac- 
cretion. As seen in the top panel, the value for (N) scales 
almost self-similarly, but with a slight increase in overall 
normalization as has been previously noted in the literature 
dGao et al.l2004l:lKlvpin et al.l2010tlGao et al.l2010t) . indicat- 
ing that more massive halos contain somewhat more subhalos, 
proportionally, than lower mass halos. 

More interesting is the trend in 0:2 with mass, shown in 
the lower panel of Figure [3] Here we see evidence that, for 
subhalos with v„ CCJll b/vmax,host ~ 0.25, self-similarity breaks 
down much more significantly, where lower-mass hosts tend 
to have systematically higher values of 012, which implies in- 
creased scatter beyond Poisson statistics. There are a number 
of possible explanations for this trend; it may be due to the 
fact that lower mass objects exist in a wider range of envi- 
ronments than higher mass halos. As shown earlier in Figure 
[5] the immediate environment surrounding a halo can have a 
significant impact on its subhalo population. High-mass halos 
always form in regions of high bias, meaning that the environ- 
mental impact on a high mass object's subhalo population is 
minimized. Low-mass objects are able to form in both bi- 
ased and unbiased regions and in regions near to higher mass 
halos, creating a larger scatter in their subhalo populations. 
While not shown, we have made analogous plots to Figure [2] 
for hosts of higher mass and confirm that these objects typi- 
cally display weaker environmental dependence in their sub- 
halo populations. 

It is also interesting to no te that previous studies, in par- 
ticular |&aytsovet]nj(|2004]) and BK10, have not seen such 
a trend. There are a number of possible explanations for 
this, which include the improved l/r'kpc force resolution and 

6 Note that, when calculating these fits, it is possible to either calculate 
(N) and a directly from the data, or to treat them as free parameters to the 
fitting algorithm. However, in practice the differences in resulting parameters 
are less than 5%. 
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FIG. 3.— Top Panel: The trend of (N) for the 

Voce distribution function 

for hosts of varying masses. This is equivalent to the subhalo mass function. 
Magenta pluses, blue diamonds, green circles, and red triangles represent 
hosts with 200, 300, 600, and 1000 km s" 1 , respectively. Bottom 

Panel: The trend of the scatter in the distribution of the number of subhalos 
above a mass threshold. The dashed line represents 012 = 1 , where the negative 
binomial distribution reduces to Poisson. 

(250/j~'Mpc) 3 volume of the Bolshoi simulation, which al- 
lows for both better statistics and more accurate tracking of 
subhalos as they orbit around their hosts. The increased vol- 
ume of Bolshoi is particularly important if the trend is driven 
by halo environment, since a larger volume is necessary to 
probe a full range of environments. Hints of this trend, how- 
ever, can be noticed in Figure 9 of BK10. While not shown, 
nearly identical trends in (N) and 010 persist if we consider 
the mass ratio v maXtSU b/v max> host, using the maximum circular 
velocity of the subhalo at z = instead of at the time of accre- 
tion. 

To summarize, in this section, we have shown that satellites 
comparable in size to the LMC and SMC are relatively rare 
in MW mass objects, occurring in roughly 9% of the poten- 
tial hosts. Their abundance is dependent on a number of host 
properties, such as the host mass and environment. The latter 
can provide a boost at the 25% level. The proximity of M31 
puts the Milky Way in a higher density environment, which in- 
creases the likelihood of an LMC/SMC pair to 12%. Finally, 
it was shown that for massive subhalos, the scatter in the ex- 
pected number of subhalos systematically decreases with in- 
creasing host mass. This may be related to the environmental 
variations previously observed; lower mass halos are able to 
form in a wider range of environments than higher mass ob- 
jects. In the next section, we turn to the task of comparing our 
simulation results with observations from the SDSS, which 
requires the adoption of an algorithm for adding magnitudes 
to simulated galaxies. 

4. SATELLITE STATISTICS FOR GALAXIES 

We turn next to a detailed comparison of these results with 
observational measurements. The classical satellite galaxies 
around the MW have been well-studied, with the SMC and 
LMC being the only satellites brighter than My = -16. In- 
deed, the next brightest object, Sagittarius, with My = -13.1, 
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is roughly 3 magnitudes dimmer than the SMC. Motivated by 
this, we search for the distribution of the number of satellites 
of My < -16 in other galaxies, both simulated and real. In 
order to distinguish the results in this section from the previ- 
ous sections, we adopt the formalism where N sats refers to the 
number of satellite galaxies around a host, identified based 
on their luminosities. This contrasts with mass-selected dark 
matter subhalos, which we characterized in the previous sec- 
tions using the nomenclature N su b s - 

4.1. The Observational Sample 

Our recent companion paper (iLiu et all |2Q10L hereafter, 
L10) studied the population of such LMC and SMC- 
luminosity objects in the SDSS. Here, MW analogs were 
identified in the SDSS main sample by locating objects with 
absolute magnitudes M r = -21 .2 ± 0.2 with no brighter com- 
panions within a line-of-sight cylinder of radius 500 kpc and 
length 28 Mpc (corresponding to a redshift of ±1000 km 
s" 1 ). 7 The abundance of satellite galaxies was then measured 
by calculating the likelihood of having N sa t s neighbors in the 
photometric sample with apparent magnitudes 2-4 dimmer in 
a fixed 150 kpc aperture. While we discuss details of their 
method below, our goal here is to compare our simulated re- 
sults to their observational constraints. 

4.2. Modeling Galaxy Luminosities 

Because internal kinematics are difficult to measure obser- 
vationally for dwarfs outside the Local Group, the observa- 
tions of L10 cannot be directly compared to the results from 
the previous section, which determined abundances based on 
the Vmax of the satellites. Instead, we employ a SubHalo 
Abundance Matching (SHAM) algorithm to our simulation 
to map magnitudes onto dark matter halos dKravtsov et al 
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2010[lBehroozi et aUteOlObTT Here, uminosities are assigned 
by matching the number density of observed objects brighter 
than a given luminosity with the number density of simulated 
objects greater than a given v max . For subhalos, v max at the 
time of accretion is used instead of the present v max (z = 0). 
We also impose a log-normal scatter between halo mass and 
luminosity. In our base model, we assume a scatter of 0.16 
dex, consiste nt with constraints on larger group and cluster 
mass objects dMore et al]|2009t iBehroozi et alJl2010bl) . 

This algorithm has been shown to be extremely successful 
in matching several statistical properties of the galaxy distri- 
bution, including the luminosity dependence and redshift scal- 
ing of the two-point gal axy correlation fun ction for galaxies 
brighter than M, < — 19 (IConrov et a l. 2006), the galaxy-mass 
correlation func tion dTasitsiom i et al. 20 04]), and galaxy three 
point statistics dMarfn et alj|2008l) . iKlvpin et al.l d2010l) have 
shown that a similar procedure can simultaneously match the 
luminosity function and the Tully-Fisher relation locally. 

Here, we match halos directly to the binned r-band lu- 
minosity fun c tion i n the local Universe as measured by 
iBlanton et all d2005l) for galaxies with M,. < -12.8. While 
there is some uncertainty in these measurements depending on 
the detailed surface brightness correction, the objects we are 
interested in here are bright enough that this effect is unimpor- 
tant. This luminosity function also reprod uces the more recen t 
number densities measured in SDSS DR7 dZehavi et al.l2010l) 

7 The value for the r-band absolute magnitud es was determined using the 
K-correction code of Blanton & Roweis 12007) on a sample population of 
500,000 SDSS galaxies. 



down to M r ~ -17. The results of this procedure are shown in 
Figure|4] Here, the red line shows the relation between v maXiacc 
and magnitude, while the blue points show v max at z = 0. The 
systematic movement to lower v max for many objects is due 
to tidal stripping of subhalos. The resulting simulated cat- 
alog has a distribution of galaxies that is complete down to 
M r = -15.3. Note that our pure-abundance matching measure- 
ments from Figure 2] gives most likely mass estimates for the 
LMC and SMC to be v max = 85 ± 6 and 65 ± 6km s" 1 respec- 
tively. However, as shown in iTrujillo-Gomez et"aTl d2010l) . 
ignoring the effects of baryons likely increases the maximum 
circular velocity by roughly 10% (see their Figure 5), bring- 
ing our SHAM mass estimates into excellent agreement with 
the observations 

Using luminosities assigned in this way, we now make mea- 
surements analogous to those of L10. We select for MW can- 
didates in the same way as L10. First, we make a mock light- 
cone by taking the z = SHAM catalog of the Bolshoi simula- 
tion, placing an observer at the center of the box, and adding 
redshift space distortions. We passively evolve galaxy mag- 
nitudes to the appropriate redshift using the S DSS constraint 
M r (z) = M r (z = 0. 1 ) + 0(0. 1 - 1), with Q = -1.3 dBlanton et alj 
120031) . We then identify all isolated objects by selecting halos 
with M, ■ = -2 1 .2 ± 0.2 (the upper grey band in Figure H] and 
excluding objects having a brighter neighbor in a cylinder of 
radius 500 kpc and length 28 Mpc. The resulting v max dis- 
tribution is shown in Figure [5] This distribution is somewhat 
lower than the best current estimates for the maximum rota- 
tional velocit y of th e MW, v ra , = 220 km s~' dXue et al.ll2008l: 
Idnedin et alJl20lOh . This difference should not be surprising 
since v max includes only the dark matter, while v ro , is the total 
rotational speed of a n actual halo that includes th e full impact 
of baryons. Further. iTruiillo-Gomez et ail (120101) have shown 
that, because the Milky Way sits near the peak of the star for- 
mation efficiency, baryons play a significant role in determin- 
ing the circular velocity profile of MW mass halos, resulting 
in a value of v max that is underestimated in dark matter-only 
simulations. However, the virial mass of the galaxies that pass 
our MW selection peaks at 1.3 x 10 12 M ^, in excellent agre e- 
ment with the estimated MW mass from lBusha et alj d2010l) . 

We select for MC-analogs using a method similar to L10, 
identifying all satellites within the virial radius of our MW 
analogs that are 2-4 magnitudes dimmer (lower grey band 
in Figure |4). The most likely v max values for our MC 
analogs are ~ 80 and ~ 70 km s " 1 , in agreement with the 
l-a constraints from observations(|van derMarel etani2002t 
IStanimirovic et al.ll2004l:lHarris & Zaritsky||2006l) . 

4.3. The Definitions of a Satellite 

Before presenting detailed comparisons between the L10 
observational data and our model, we note that the definition 
of a "satellite" differs between the two analyses. Because L10 
used a photometric sample with hosts selected on luminosity 
and no spectroscopy for their satellites, the number of satel- 
lites was defined as the number of photometric pairs around 
a host. Specifically, they identified isolated MW-luminosity 
galaxies in the SDSS spectroscopic sample and counted the 
number of photometric pairs with Am r = 2-4 within a fixed 
angular separation, nominally taken to be 150 kpc. They then 
applied a background subtraction to account for projection ef- 
fects and correlated structures. While our SHAM catalog will 
let us select isolated host objects in an analogous way, be- 
cause modeling the full distribution of galaxy apparrent mag- 
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FIG. 4. — The relation between magnitude and v max for simulated ha- 
los+SHAM. The red line represents the relation between M r and v acc , the 
maximum circular velocity at the time of accretion, which was used to se- 
lect the magnitude. The black curve and blue points show the relation the 
between magnitude and Vmax at z — (as opposed to Voce)' The upper grey 
band shows our magnitude-criteria for selecting Milky Way-like hosts, while 
the lower gray band represents our range for selecting Magellanic Cloud-like 
satellites. 
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FIG. 5. — The distribution of Vmax f° r simulated galaxies, as selected by 
luminosity. The black histogram represents galaxies that pass our Milky Way 
selection, with M r = —21. 2 ±0.2, while the blue and red histograms represent 
Vmax and v a cc for our Magellanic Cloud-like subhalos that are 2-4 magnitudes 
dimmer than their hosts. The dotted line shows the preferred value for v,.„ of 
the MW. 



nitudes is beyond the scope of our model, we do not perform 
an identical calculation. Instead, we can select an appropri- 
ate definition of a satellite galaxy to mimic the observational 
effects. 

L10 presented two methods for background subtraction, 
which were used to calculate the number of satellites around 
a host galaxy. The first, which we refer to as isotropic sub- 
traction, assumed that background objects had, on average, 
no dependence on position (implicitly ignoring the impact of 
correlated structures), to calculate an average background for 
their entire host population. Using this assumption, it is possi- 
ble to calculate the average background to very high statistical 
accuracy. Their second method, which we refer to as annulus 
subtraction, measured the background object counts around 
each galaxy by considering the region just outside the pro- 
jected region surrounding the galaxy. For completeness, we 



FIG. 6. — The 3D cross-correlation function, £( r )MW-MC, between Milky 
Way like hosts and all LMC and SMC magnitude satellite galaxies. The black 
dashed line represents a power-law fit to £(r). For reference, the dotted red 
line denotes the auto-correlation b etween galaxies with — 2 1 < M r — 5 log(/i) = 
-20 in SDSS DR7 as measured in Zehavi et al. 1 2010). 

make comparisons to both methods. 

Because it removes the effects of correlated structures, 
comparison to the annulus subtraction method of L10 is 
straightforward. We define as a satellite galaxy any object 
within the appropriate spherical volume (taken to be 150 kpc 
in L10) that has an absolute magnitude that is 2-4 mag dimmer 
than its host. To compare to the isotropic subtraction method, 
which neglects the impact of correlated structures along the 
line of sight (see Section 3.2 of L10 for further discussion) 
requires additional modeling. We account for this by defining 
a satellite to be any object inside a cylinder of radius 150 kpc 
and length ro, where ro is the correlation scale for the Milky 
Way host-satellite analog. 

This cross-correlation function between MW-luminous 
galaxies and objects 2-4 magnitudes dimmer is plotted in fig- 
ure^ and is well fit by a power law, 



ttr) = (r/r r 



(5) 



with ro = 5.3 Mpc and 7 = 1.7, with a similar correlation 
length to measure ments for all galaxie s of MW-like magni- 
tudes inferred from lZehavi et al.l ( 120101) . These authors mea- 
sured ro = 5.98 Mpc//i and 7 = 1.9 for the autocorrelation 
function of galaxies with M r = -20.7 to -21.7. We there- 
fore consider a cylinder with radius 150 kpc and length 5.3 
Mpc (the correlation length) to be the appropriate aperture for 
satellite selection to compare to the isotropic background sub- 
traction of L10. 

Coincidentally, this correlation scale is very s imilar to the 
size of the expected redshift-space distortions. lEvrard et alj 
d2008l) showed that a ~ 10 12 M Q halo has a typical velocity 
dispersion of 125 km s" 1 ; at z = 0.08, where the bulk of our 
simulated halos live, this corresponds to a 3.7 Mpc spread due 
to the redshifts of the satellite population. A cylinder of this 
length thus represents the most accurate measurements one 
could make for the number of satellites around a host using 
spectroscopic data (without including some sort of additional 
background subtraction). 

In Figure [7] we compare the P(N su i, s ) distribution for four 
different satellite definitions: spherical and cylindrical aper- 
tures (which correspond to the L10 results using their annulus 
and isotropic background subtractions methods, respectively) 
as well as selecting satellites based on v max and abundance- 
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FIG. 7. — The probability of having N sats LMC and SMC-type objects in 
simulations based on SHAM and v max for various satellite definitions. Red 
circles count objects with v max > 50 km s~' within the virial radius; the black 
diamonds, blue triangles, and green boxes count objects 2-4 magnitudes dim- 
mer than their host using a SHAM model with 0.16 dex scatter. The black 
diamonds were calculated using all objects inside the virial radius. The blue 
triangles and green boxes are equivalent to the isotropic and annulus back- 
ground subtractions used in L10: a cylindrical aperture of radius 150 kpc and 
length 5.3 Mpc (blue triangles), and a fixed 150 kpc spherical aperture (green 
boxes). 

matched magnitude. Here, the red circles consider subhalos 
within the virial radius of their host that have v max = 50- 80km 
s -1 . The black diamonds, blue triangles, and green boxes all 
use magnitude-selected satellites, counting objects 2-4 magni- 
tudes dimmer than their host using the aforementioned SHAM 
model with cuts analogous to the observational criterion of 
L10. The red diamonds represent the "correct" number of 
satellites, using all objects inside the virial radius, The blue 
triangles use a cylindrical aperture, analogous to the isotropic 
background subtraction of L10, while the green boxes use 
a spherical aperture, equivalent to the annulus subtraction 
method of L10. All error bars were calculated using the jack- 
knife method with 10 sub-regions. 

A number of trends are readily apparent for the relations 
in Figure [7] First, there is generally good agreement between 
the P(N su bs) distribution for objects based on v max (red circles) 
with our SHAM magnitudes (black diamonds). This is to be 
expected if SHAM preserves the observed mass-luminosity 
relation in this regime, something that will only happen if an 
appropriate mass function is used. It is also interesting to note 
that our calculations for the cylindrical aperture (blue trian- 
gles) always lie above the other measurements for N su b s > 1 . 
This shows that correlated structures are an important effect 
which is typically comparable to the signal from satellites that 
are formally inside the halo's virial radius. Our measurements 
using a fixed 0.15 Mpc cylinder (green box) lie below other 
measurements for higher values of N sats . This should be un- 
surprising, since the typical virial radius for our selected halos 
is roughly 250 kpc. It is interesting to note, however, that the 
cumulative probability for hosting (at least) one satellite dif- 
fers by less than 5% when calculated using the virial radius 
instead of a radius of 150 kpc, while this difference increases 
rapidly for higher values of N sats . All of these trends indicate 
the importance of careful background subtraction to eliminate 
this signal (such as the annulus subtraction of L10), or other- 
wise of comparing directly to simulations with similar selec- 
tion criteria so as to properly interpret results. 

As concerns the length of our cylindrical aperture, we ob- 



tain almost identical results for any length beyond ~ IMpc, 
beyond which the correlation function falls off substantially. 
Varying the cylinder length from 2-8 Mpc changes our results 
by roughly five percent, comparable to the size of our statis- 
tical errors. Thus, the comparison between satellites defined 
using our cylindrical aperture and the isotropic background 
subtraction technique of L10 should be robust. 

The numerical values for many of these probabilities — se- 
lection by r V i r , as well as our cylindrical and fixed spherical 
apertures — are presented in Table [2] Note that, while we 
are only counting objects as satellites if they are 2-4 magni- 
tudes dimmer than their hosts for the most accurate compari- 
son with L10, these numbers are essentially identical to con- 
sidering a threshold sample of objects 0-4 magnitudes dim- 
mer than their (MW-like) hosts. Switching to this definition 
creates only percent-level changes, well within the statistical 
error bars. 

4.4. Comparisons Between Simulations and Observations 

Figure [8] and Table [2] show the main result of our analy- 
sis, making direct comparison with the results of L10. In the 
top panel, we compare the L10 P(N sats ) measurements using 
their isotropic background subtraction (black asterisks) with 
our SHAM results using subhalos inside a fixed cylinder with 
radius 150 kpc and length 5.3 Mpc (blue triangles; cylinder 
column in Table 2) that closely mimics their selection crite- 
ria. The bottom panel shows the comparison between the L10 
annulus subtraction (black and gray asterisks) and our mea- 
surement of Ngats within a 150 kpc spherical aperture (green 
squares; Sphere column in Table 2). The black and gray aster- 
isks in the lower panel represent the range of measurements 
based on uncertainty in selecting an optimal background sub- 
traction. The difference between these sets of points thus 
gives an estimate for the size of the systematic errors in the 
measurements. Note that, after making their background sub- 
traction, L10 also calculate a systematic adjustment due to 
catastrophic photo-z failures (see section 4.2 of L10). While 
we have included those corrections in Figure [8] we have ex- 
plicitly kept the lower error bar consistent with zero for any 
measurement that was consistent with zero at the 2-er level 
before this systematic correction was included. For reference, 
Table |2] also gives the P(N sa ts) distribution for proper subha- 
los, i.e., objects within r„> that are 2-4 magnitudes dimmer 
than their host. 

The agreement between simulations and observations in 
Figure [8] is generally very good. We take this to indicate a 
success for both the ACDM and SHAM models, particularly 
indicating, with excellent statistics, that the HOD predicted 
by Bolshoi for Milky Way-like objects at the massive end 
matches observations, and there is no evidence for either a 
missing or excess satellite problem for these objects. Consid- 
ering our comparison between the cylindrical/isotropic sam- 
ple, at no N sa t s do the two results differ by more than ~ 2.5a, 
and in most situations agree to better than 1 a. The com- 
parison fares slightly worse for the spherical/annulus sample. 
Here, the simulations underpredicts the number of N sa ts = 
systems, while over-predicting the N sa ts = 1 systems at the 3 
a level. While this could indicate a real disagreement of the 
simulation and galaxy formation model with the data at this 
mass scale, there are a number of other possibilities. First, it 
is possible that the annulus-subtraction method of L10 over- 
corrects the background subtraction. This has been tested 
in simulations, using their algorithm to estimate the number 
of objects within 150kpc around our halos centers using our 
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FIG. 8. — Comparisons between the P(N sa t s ) distribution of SMC-luminous 
satellites around MW-luminous hosts in the Bolshoi simulation and the obser- 
vational results of L10. Top Panel: Comparison of the isotropic background 
subtraction model from L10. Black asterisks represent the abundances from 
SDSS hosts. The blue triangles denote the abundance of objects inside an 
aperture that most closely corresponds to the isotropic subtraction model: a 
fixed cylindrical aperture with radius 150 kpc and length 5.2 Mpc. Lower 
Panel: Same as the upper panel, but now comparing the annulus subtraction 
model from L10 (black and gray asterisks) to the simulation-equivalent of a 
fixed 150 kpc spherical aperture for identifying satellite galaxies around their 
hosts (green squares). For the L10 data, the black points represent our fidu- 
cial measurements, while the gray points add a systematic uncertainty to the 
background subtraction. 



cylindrical counts. The result, is an over preduction of the 
Nsais = term by 6% and an underprediction of the N sats = 1 
term by 5%, larger than the estimated systematical uncertain- 
ties in Table 1 of L10 (the changes in the higher N sals terms 
are of similar size to their estimated systematic uncertainty, 
and are therefore well characterized by the gray points in Fig- 
ure [SJ. Such corrections would bring the these terms within 
2- and 3- a of our simulated estimates. It is also important 
to note that the predictions from simulations depend on the 
assumed scatter between mass and luminosity. This will be 
further discussed below. 

Looking at the SDSS/Bolshoi comparisons in more detail, 
we can further investigate the dependence of the P(N sats ) dis- 
tribution on the aperture size, as done in L10. For this, we 
determine the probability of finding N sats satellite galaxies in 
a spherical aperture around their host, where the radius is var- 
ied from 100 to 250 kpc using the annulus background sub- 
traction. The results are shown in Figure [9] Here, the red 
circles, orange squares, magenta triangles, and cyan inverse- 
triangles represent the probability of having 0, 1, 2, or 3 sub- 
halos. Open symbols represent SDSS measurements, while 



TABLE 2 

Probability of hosting N sat s satellite galaxies brighter than 

Mrjmst +4 AROUND A MlLKY WAY MAGNITUDE HOST HALO WITHIN 
VARIOUS APERTURES FOR BOTH OUR SIMULATED RESULTS AND THE 
OBSERVATIONAL MEASUREMENTS FROM L 1 0. 





r a - 
vir 


Cylinder 6 Isotropic 


Sphere'' Annulus 1 " 


P(0) 

P(l) 

P(2) 
P(3) 
P(4) 


55 ±4% 
28 ± 3% 
11 ±1.4% 
3.8 ±0.7% 
1.55 ±0.4% 


55 ±5% 661}| 
28 ±3% 22+[| 
11 ±2% 
3.8 ±1.2% 1.5^ : ° 

1.5 ±1.1% 0.9+°<: 7 


71 ±3% 8l!|j 

23 ±2% 
5.2 ±1.6% 3.5+JJ 
1.0 ±0.8% 1.6^? : ? 
0.2 ±0.1% l.l^f 



"halos contained within the virial radius of the host 

*halos within a cylinder with radius 150 kpc and length 5.2 Mpc 

C L10 results counting objects within a projected 150 kpc and an isotropic 

background subtraction 

''halos within a sphere of radius 150 kpc 

C L10 results counting objects within a projected 150 kpc using an annulus 
background subtraction 
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FIG. 9. — Abundances of LMC and -type objects in the simulation, as 
a function of search radius. In all cases the filled symbols were calculated us- 
ing Bolshoi+SHAM, while the open symbols represent SDSS measurements 
from L10. Here, the red circles, orange squares, magenta triangles, and cyan 
inverse-triangles represent the probability for a Milky Way like galaxy to have 
0, 1, 2, or 3 subhalos, respectively. For the Bolshoi objects, the probabilities 
were calculated using a spherical aperture with varying radius. The SDSS 
points used the annulus subtraction method of L10. The hatched lines repre- 
sent the probabilities for hosting N sa ts subhalos inside the virial radius of the 
host. 



filled symbols come from our Bolshoi+SHAM model, while 
the hatched lines represent the P(N SMS ) distribution for satel- 
lites within the virial radius of their hosts. Again, we see good 
agreement between simulations and observations, although 
we note the persistence of higher values in the simulation for 
the N sats = 1 measurement. However, the other values of N sats 
are in much better agreement. 

One of the appeals of the SHAM algorithm is that it has so 
few free parameters, in this particular case just the scatter (al- 
though there are additional implicit assumptions, e.g., about 
subhalo stripping and star formation after satellite accretion). 
Figure[10]shows how the scatter impacts this result. Here, we 
present the probability distribution for an isolated host with 
M r = -21.2 to host Nsats satellites brighter than -16.5 as we 
vary the scatter in the SHAM model from to 0.3 dex using 
our cylindrical method for selecting satellites. While none of 
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FIG. 10. — Dependence of satellite probabilities on scatter in the rela- 
tionship between halo mass and galaxy luminosity. Here, we compare the 
P(N s ats) distribution calculted using our cylindrical method to the L10 mea- 
surements using an isotropic background subtraction. In each case, the model 
assumes a constant value of a corresponding to the scatter in log L at a given 
Vacc. Larger scatter values tend to have more systems with many satellites, 
largely because more high mass hosts are scattered into the sample. 

these models are ruled out by the data (L10 isotropic model), 
a number of systematic trends are present. In particular, there 
is a systematic increase in the probability of hosting 3 or more 
satellites, along with a decrease in the probability of hosting 

1 satellite (by up to 5%). This trend is caused not by scatter 
in the mass-luminosity relation of the satellite galaxies, but 
in the mass-luminosity relation of the host objects. As scat- 
ter increases, higher and higher mass hosts scatter into our 
magnitude-selected sample, bringing with them their richer 
subhalo populations. As we shall discuss in the next section, 
this has the effect of creating a tail to the negative-binomial 
distribution that sets the probability for a host to have N sars 
satellites. 

It is also worth noting that the dark blue points represent- 
ing a scatter of a = 0.3 dex look to be violating the trend of 
an increasing high-N sals probability. This is not due to any 
special behavior in this regime, but more represents the lim- 
its of our simulation resolution. As scatter increases, objects 
with M, ■ -5log/i > -16.5, which nominally have v acc = 90 
km s" 1 using a scatter-free SHAM, have infall masses as low 
as 65 km s"'in our a = 0.3 dex model. Because tidal strip- 
pi ng can reduce the v ma y of a subhalo by a factor close to 

2 dWetzel & Whitdl201oT iGuo et al.lloTot) . this puts such ob- 
jects dangerously close to the Bolshoi completeness limit of 
50 km s" 1 . It should be noted, however, that the completeness 
of subhalos in the Bolshoi simulation is highly dependent on 
the host mass. Because of stripping of subhalos after accre- 
tion, we expect to be losing potentially 10-15% of the satellite 
population around ~ 10 12 /T 1 M Q mass halos with scatter this 
high. It's important to appreciate the high resolution neces- 
sary to fully model a galaxy population using SHAM. While 
a scatter of a = 0.3 is likely ruled out for massive galaxies 
dBehroozi et al.ll2010bl) . if such scatter were appropriate the 
Bolshoi simulation would only be able to populate galaxies 
down to M r ~ -17.5 while losing fewer than 10% of the ob- 
jects. 

Finally, we should note that there is no compelling reason 
to assume that both massive and low-mass galaxies exhibit the 
same scatter in the mass-luminosity relation. We have inves- 
tigated this by quantifying the variations in the P(N sa ts) distri- 



bution where we keep the mass-luminosity scatter of the hosts 
fixed at 0.16 and vary the scatter of the subhalos from 0-0.3. 
While we see similar trends with, e.g., the N sars = 1 term drop- 
ping with increased scatter, the effect is significantly weaker 
than shown in Figure [10] where we allow the scatter in the 
host mass-luminosity relation to vary. While it is very pos- 
sible that the scatter should be even higher than 0.3 for our 
satellite galaxies, the 50 km s" 1 resolution limit of Bolshoi pre- 
vents us from exploring a wider range. Thus, over the ranger 
we are able to probe, the scatter in the host mass-luminosity 
relation, which sets the range of host masses sampled, has a 
larger impact that the scatter in the satellite mass-luminosity 
relation. 

4.5. Environmental Effects 

We now revisit the discussion rega rding the environmental 
impact on the subhalo population of 33.21 using observation- 
ally motivated environmental measures. Recall that in 33.21 
we defined the local environment as the dark matter density 
in collapsed objects larger than 2 x 1O H M0. Here, we pre- 
form a similar analysis using luminosity as an environmental 
proxy. Specifically, using the isolated central galaxies from 
our SHAM catalog as selected above, we measure the total lu- 
minosity from galaxies brighter than M r = -20.8 (which corre- 
sponds to M r - 5 log(/i) = -20 in the usual SDSS units) within 
a sphere of 1 Mpc around our hosts, which we take to be our 
environment proxy. We then split the sample into our most 
and least luminous (i.e., dense) quartiles and re-measure the 
P(N sa ts) distribution for satellites inside their host virial radii. 
The results are shown in Figure [TT] Here, the black points 
represent the distribution for all isolated hosts, while the red 
and blue represent those in high/low density environments. 
As before, we see a clear offset for galaxies in different envi- 
ronments. While the differences are only at the 1-2 sigma lev- 
els, high-density halos are systematically more likely to have 
massive subhalos than those in low density environments. In 
particular, the underdense halos are 14% more likely to host 
no massive subhalos than their dense counterparts. Similarly, 
halos in dense environments are twice as likely to host two or 
more subhalos than those in low density environments. This 
effect should therefore be an observable manifestation of as- 
sembly bias. 

5 . STATISTICS OF SATELLITES IN LUMINOSITY-SELECTED HOSTS 

As noted previously, the inclusion of scatter in the halo 
v max -galaxy luminosity SHAM creates a significant tail in the 
probability distribution for a host to have N snts satellite galax- 
ies. Consequently, the negative binomial distribution that was 
found to describe the P(N su t> s ) distribution for the subhalos in 
a mass-selected host does not produce a good fit to P (iV sa ts) 
when applied to luminosity-selected galaxies. Additionally, 
the change in star formation efficiency with host halo mass 
breaks the self-similar ity that was observed for host halos 
dBehroozi et alJl2010bl and references therein). In this sec- 
tion, we quantify the P{N sats ) distribution for satellite galaxies 
around a luminosity-selected host. 

In Figure [12J we show a series of probability distribu- 
tions for a M r = -20.5 galaxy to host N sats subhalos brighter 
than -19.8, -18.3 and -16.8 (equivalent to M r -51og(/z) = 
-19,-17.5,-16). This plot was generated using a scatter of 
a = 0.16 dex. A low probability tail clearly extends to high 
N S ats- Such a tail was not present w hen c onsidering v max - 
selected hosts and subhalos in section [331 This is empha- 
sized by the dotted lines, which present the negative binomial 
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FIG. 1 1 . — The probability distribution for a host galaxy to host NS satel- 
lites within its virial radius split by environment. The black diamonds rep- 
resent the distribution for all isolated hosts (and are identical to the black 
diamonds in Figure|7J. The red and blue triangles represent the distribution 
for hosts in high and low density environments. Here, we use the total lumi- 
nosity in galaxies brighter than M r = -20.8 as our proxy for density. 




5 10 15 



FIG. 12. — Probability of hosting N sa ts satellite galaxies in a Milky Way- 
luminosity host, for various satellite selection criteria. The red, green, and 
blue points represent satellites brighter than M r = — 19.8. — 18.2, and —16.8. 
Dotted lines represent the best fit negative binomial distribution, while the 
dashed line shows the best fit negative binomial plus power-law distribution. 

distributions, equation|3] that characterized the P(N su bs) distri- 
bution. While the negative binomial provided an excellent fit 
to the data for our mass-selected objects, it clearly fails here. 
While this tail may not appear to be a serious concern, since 
it is not important until P(N sats )~ 1%, large surveys such as 
the SDSS contain more than 10,000 massive clusters, so even 
these small effects may be visible. 

As a better fit, we present a modified binomial distribution 
by adding a power law tail, 

P(N\r,p,T,T)= r(jV+r) p r (l- p f+Te- TN . (6) 

Here, the first term is just the negative binomial, with r and 
p functions of 0:2 and (N), as defined in Equation [4] The 
final term models the high-A^,,, tail that comes about as a re- 
sult of the v max -luminosity scatter, where T and r control its 
slope and amplitude. The dashed lines in Figure Q~2] show the 
best fits to Equation [6] This clearly provides a significantly 



improved fit to the simulation measurements, albeit at the ex- 
pense of adding two new fit parameters. Note, however, that 
the fit tends to have problems in the cross-over region and that 
the transition to the high-N sats tail is somewhat more gradual 
than this model is able to reproduce. While we leave this for 
future investigations, it is worth noting that changing the scat- 
ter in our SHAM model has a direct impact on T and r. High 
values of scatter will have a more pronounced tail, which will 
cause both the amplitude and slope of this tail to increase. For 
simplicity, for the rest of this work, we will limit ourselves to 
using a model with a = 0.16 dex. 

While we could present our best-fitting parameters to the 
distributions in Figure |T2] a much more useful process is to 
model the P(N sals ) distribution as a function of both host mag- 
nitude and limiting satellite mag nitude, as has bee n done for 
mass-selected subhalos (see, i.e., iGao et all 12004ft . Unfortu- 
nately, this proves to be much more difficult for galaxies than 
for halos. Because the dark matter power spectrum is mostly 
featureless in the relevant regime, large halos are nearly self- 
similar versions of smaller halos. As an example of this, the 
top panel of Figure [3] shows the mean number of satellites 
above a given mass cut for mass-selected host halos. This dis- 
tribution exhibits a power-law behavior with nearly identical 
parameters regardless of the host mass. This is not the case 
for magnitude-selected objects. 

Figure Qj] presents the best-fit parameters for the satel- 
lite P(N sats ) distribution as modeled by Equation [6] around 
galaxies as a function of limiting subhalo luminosity (i.e., 
magnitude difference). The red, orange, yellow, dark 
green, green, and blue points represent hosts with M r = 
-20.3,-20.8,-21.3,-21.8,-22.3, and -22.8. As can be 
clearly seen in the top left panel, which shows the trend of 
the mean number of satellite galaxies ((N)) with luminosity 
ratio, galaxies do not scale self-similarly. At a given lumi- 
nosity ratio, brighter galaxies are significantly more likely to 
have a larger N sats (expected number of satellites) than dim- 
mer galaxies. 

This is a manifestation of the varying star formation effi- 
ciency with halo mass as seen in, e.g., Figure [4] For masses 
above 10 12 M Q , star formation efficiency becomes less effi- 
cient, so the rate of increase of M r with v max is shallow. At 
masses lower than 10 12 M , star formation efficiency drops 
rapidly with decreasing mass, leading to a much steeper rate 
of change of M r with respect to v max . Thus, for a fixed bin in 
M r , the corresponding logarithmic range in v max will be much 
less for low-mass galaxies as compared to high-mass galaxies. 
As such, while the (N) values for mass-selected subhalos lie 
on a single line in Figure [3] luminosity-selected halos do not 
scale in such a manner, resulting in the different lines in the 
top-left panel of Figure Q~3] brighter host galaxies have more 
satellites at a fixed luminosity ratio, L sat /Lh os t, than dimmer 
hosts. 

A number of the trends present in Figure [13] can be readily 
understood in terms of simple models. As discussed above, 
the behavior of (N) is simply a manifestation of the chang- 
ing star formation efficiency with halo mass. Similarly, the 
trends in r can be understood though simple scaling relations. 
First, consider a fixed luminosity ratio, \M r j, 0St -M r ji m \. Since 
brighter hosts correspond to more massive halos, these ob- 
jects live in a steeper part of the mass function, where a fixed 
scatter creates a broader selection of host masses, the larger 
of which are more likely to host bright satellites. Because of 
this, the exponential tail falls off more slowly, resulting in a 



12 



BUSHA ET AL. 



A 

£ 4 
v 



■ i ■ ■ - i i i ■ 



■i i ■ 



\ 

■■ ^ ~ at sfc**=*>.s 



-6 -5 -4 -3 -2 -1 

M r w-M rJim 





-6 -5 -4 -3 -2 -1 



-6 -5 -4 -3 -2 -1 



1.00 r 




-6 -5 -4 -3 -2 -1 

M r , host - M rJim 



FIG. 13. — The fit parameters specifying the P(N) distribution (equation [6) as a function limiting satellite magnitude, represented as a magnitude difference 
(luminosity ratio) Different colored points represent different host populations: violet, blue, green, dark green, orange, and red points are hosts with M r = 
-19.8,— 20.3,-20.8,— 21.3,— 21.8, and —22.3. Dashed lines represent out fits to these distributions, as given in Equation[7J 



lower t. The trend of r increasing with luminosity ratio for a 
fixed host mass. 

Guided by these observations, and noting that a.2 and T 
don't have strong trends with either host magnitude or limit- 
ing luminosity ratio, we propose a 7-parameter fit to generally 
describe the full P(N\M r i JOS ,,M r n m ). We use our generalized 
negative-binomial plus power law as defined in Equation |6j 
and keep an, T, and r constant. However, we parameter- 
ize (N) as a linear function, the slope of which depends on 
M r fost-M r u m and the offset of which depends on host mag- 
nitude, M r host. Our full fit for Equation [6] then becomes a 
function of the parameters NQ. a ,N^b,N\ ia ,N\fi,a2,TQ, and to 
defined by: 



P(N\M r j wst ,M r j im ) = 



where 



T(N+r) 
r(r)r(JV+l) 



/(I 



-pf+T e 



toN 



(7) 



1 



1 

a? — 1 



l){N}(M rMos ,,M r Mm) 



(8) 



(N)(M r ,host,M ri i im ) 



: _ J QNo.fl+No.b *log(-M rMsl +5 log(A)) _ 

(M rto -M Wim )10 Wl ^ +iVl ' i '* log( - M * ! ' +51og(A)) . 



We fit this functional form to the measured P(N sa i s ) distribu- 
tion for all host halos down to M r = -18.5 and their satellites 



TABLE 3 

Parameters specifying P(N) distribution in Eq.[8] 



parameter best-fit value 



N .a 

N , b 
Ni, tt 

OL2 
T 



-53.70 
40.11 

-39.05 
29.25 
1.3 

1.3 x 10-* 
0.023 



down to M r = -16 in Bolshoi stacked in bins of 0.5 mag. We 
present our best fitting parameters for Equation [8] in table [5] 
As a final validation of this result, Figure [14] shows the mea- 
sured P(N) distributions for a range of M r j ms , and M r ji m , along 
with the fits from our model in Eq. [7] Note that, because the 
slope and intercept of (N)(M rt hast >M r ji m ) are linear functions 
in log space, the exact value of (N) is highly sensitive to the 
parameters No. a ,No.b,N\. a , and N\ t b- Because of this, it is nec- 
essary to constrain these parameters to 4 significant figures 
(which is easily done given the constraints from the Bolshoi 
simulation). 

While our nominal fit can be improved using a r that is a 
power-law in M r j, osl —M r ,iim with a normalization that depends 
on host mass, the actual fit does not significantly improve the 
reduced x 2 , so we elect to hold it constant. 
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FIG. 14. — Probability distributions for galaxies to host N s ats satellites (diamonds). The four panels represent different selections based on host magnitudes 
M r ,host - -19.5,-20.5,-21.5,-22.5 for the top left through lower right. The blue, green, and red points represent various limiting satellite magnitudes, M, ■;,„, = 
-16.5,-17.5,-18.5. The dashed lines represent our model with fit parameters from Eq.|7]and Table[5] 



6. CONCLUSIONS AND SUMMARY 

In this work, we have examined the likelihood for simulated 
MW-like dark matter halos to host substructures similar to the 
Magellanic Clouds. We have performed this analysis using 
both dark matter kinematic information and an abundance- 
matching technique which assigns galaxy luminosities to re- 
solved dark matter halos and allows direct comparison with 
observations (e.g., L10). The main result of this work is that 
MW-like objects have a 5-11% chance to host two subhalos 
as large or as luminous as the SMC, in basic agreement with 
previous simulation results (e.g., BK10), as well as with ob- 
servations. While we have extensively compared our results 
for the full probability distribution for a MW-luminous galaxy 
to host N sars satellite galaxies to the observational measure- 
ments from LI 0, finding good agreement. This results are also 
in qualit ative agreement with previous results u sing smaller 
samples (IChen et al.ll2006tlJames & Tvoryll2010l) . 

Our results have a number of implications. First, this is 
a validation of the ACDM paradigm down to the level of 
10 10 /2 -1 Mo objects. At this level, there is no indication that 
either the CDM paradigm or our standard picture of galaxy 
formation, including the relatively tight relationship between 
galaxy luminosities and halo masses, breaks down. In partic- 
ular, it appears that the statistics of satellite galaxies at this 



mass in LCDM are in very good agreement with observations 
which sample a wide range of galaxy environments. 

One possible criticism of the model we have used for con- 
necting galaxies to halos (SHAM) is that it may be so robust 
as to always reproduce the observed statistics, given its as- 
sumption of the correct global luminosity function. In order to 
evaluate this, it is important to keep track of the exact assump- 
tions of the approach and how these relate to the particular 
statistics measured. First and foremost, SHAM assumes that 
iWx of a halo or subhalo (at accretion into a larger system) 
is the only parameter that sets the luminosity of the galaxy 
it hosts. While this is robust to, e.g., a temporally-variable 
star formation efficiency with halo mass, it would differ from 
model in which there was environmental dependence to the 
galaxy formation recipe in addition to the inherent environ- 
mental dependence of the halo mass function. For example, 
if subhalos of a given mass evolve very differently than halos 
of the same mass in the field, our model would not produce 
good agreement. Additionally, if our simulated dark matter 
halo mass function was significantly off from reality (due to, 
for example, an incorrect a% or f2 m ), the masses of the MCs 
as estimated from abundance matching would not provide a 
good match to the direct kinematic measurements, as we find 
here. Finally, this is a test of the SHAM treatment of subhalos 
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in a mass regime lower than has been previously explored. 

We do see some manifestations of assembly bias with re- 
gard to the dark matter subhalo population inside MW-like 
galaxies: halos in higher density regions host more massive 
subhalos than those in lower density regions. This effect ap- 
pears to be most pronounced for relatively local measure- 
ments of density, and suggests that the MW's proximity to 
M3 1 boosts the likelihood for us to see a LMC/SMC pair by 
about 25%. The effect is strongest for densities measured 
on the scale of 1 Mpc. We note that such a boost seems to 
be present if we use either the total local mass or the total 
local luminosity as an environmental measure. This should 
make the effect an observable manifestation of assembly bias 
(though detecting it will require careful treatment of corre- 
lated structure). It is interestin g that this result is in qu alita- 
tive agreement with the work of llshiyama et al.l (l2009bl) . who 
saw a similar trend with environment but claimed it was most 
sensitive at the ~ 5/i _1 Mpc scale. They note that we are in an 
underdense region on these scales, and cite this as a possible 
ingredient for solving the missing satellite problem. While 
our simulations do not allow us to determine the probabili- 
ties of finding objects smaller than the Magellanic Clouds, if 
our observed trend were extrapolated down to lower masses it 
would result in an opposite effect, where the presence of M3 1 
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could serve to exacerbate the missing satellite problem. 

Finally, we have generalized our measurements for the 
P{N sats ) distribution over a wide range of host and satellite lu- 
minosities. Our results show a significant deviation from the 
Poisson expectation for larger values of N sa t s - While these de- 
viations are typically present only in the low-likelihood tail of 
the distribution, the large volume of surveys like SDSS should 
allow tests of these predictions with observations. Deviations 
from Poisson statistics in the number of satellites may also 
have consequences for HOD modeling, which currently as- 
sume a Poisson distribution for the number of galaxies given 
the halo mass. 
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